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ABSTRACT 



We introduce a class of fourth order symplectic algorithms that are ideal 

^ for doing long time integration of gravitational few-body problems. These algo- 

' rithms have only positive time steps, but require computing the force gradient 

^ ■ in additional to the force. We demonstrate the efficiency of these Forward Sym- 

Tj- ■ plectic Integrators by solving the circularly restricted three-body problem in the 

^ ' space-fixed frame where the force on the third body is explicitly time- dependent. 

O ■ These algorithms can achieve accuracy of Runge-Kutta, backward time step and 

■ corrector symplectic algorithms at step sizes five to ten times as large. 

O 

?-j I Subject headings: symplectic integrators, higher order, trajectory computation, 

^ I long time simulation 

> 

X 

H ' 1. Introduction 

Symplectic Integrators (SI) (Yoshida 1993; Channell and Neri 1996) conserve all Poincare 
invariants when integrating classical trajectories. For periodic orbits, their energy errors are 
bounded and periodic, in contrast to Runge-Katta type algorithms whose energy error in- 
creases linearly with integration time (Gladman, Duncan and Candy 1991). Symplectic 
algorithms are therefore better long time integrators of astrophysical many-body problems 
(Wisdom and Holman 1991; Saha and Tremaine 1992) and have been studied extensively 
in the literature (Candy and Rozmus 1991; Yoshida 1993; McLachlan 1995; Koseleff 1996; 
Blanes et al. 1999). However, current higher order SI seem to suffer two disadvantages. First, 
the number of force evaluations require by these algorithms proliferates much more rapidly 
than non-symplectic algorithms. For example, a six order Runge-Kutta-Nystrom algorithm 
(Albrecht 1955) requires only five force evaluations; a six order symplectic algorithm with 
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negative intermediate time steps requires at least seven force evaluations (Yoshida 1990). 
Second, they seem to have a shorter stability range. This may also due to the use of nega- 
tive time steps for cancelling lower order errors. In this work, we consider a class of fourth 
order Forward Symplectic Integrators (FSI), which have only positive time steps (Chin 1997; 
Chin and Chen 2002). (We have referred to these as gradient symplectic algorithms previ- 
ously.) The price for having only positive time steps is that one must compute the force 
gradient in additional to the force. This additional effort for gravitational few-body problems 
is very modest, but the resulting gain in numerical efficiency and stability is tremendous. 
We will show below that these FSI can achieve accuracy of existing symplectic integrators 
at step sizes five to ten times as large. These very powerful algorithms should be considered 
by coUeauges doing long time few-body simulations. 

Despite the fact that Wisdom and Holman (1991) have questioned in 1991 why there 
were no FSI beyond second order, this class of algorithms has been unrecognized and over- 
looked until recently. The reason for this is simple; classical dynamics is time reversible, there 
is no compelling demand for purely positive time steps. The demand for positive time steps 
only arises when one has to solve time-irreversible equations, such the diffusion (Auer et al. 
2001), Fokker-Planck (Forbert and Chin 2001), and Navier-Stoke equations. However, once 
FSI have been derived, they can be used to solve both time-reversible and time-irreversible 
evolution equations. For time-irreversible systems, they are the only decomposition algo- 
rithm possible, but even for time-reversible equations, they give rise to algorithms of great 
efficiency and stability. 

In the next Section, we will describe in detail how this class of algorithms fits into the 
large picture of symplectic integrator development. We will show that our FSI traces it 
lineage to Ruth's seldom used third order algorithm (Ruth 1983). In Section 3, we describe 
our fourth order FSI and the one-parameter family of integrators. In Section 4, we compare 
various algorithms by computing a lengthy and intricate closed orbit of a circularly restricted 
three-body problem. Physical three-body problems, such as the Sun-Earth-Moon, or Sun- 
Jupiter-Saturn, are too "tame" for discriminating the merit of these very powerful algorithms. 
Here, we also note a puzzling fact that some algorithms can appear to converge at higher 
order than expected. We summarize our conclusions and suggest some future direction of 
research in Section 5. 
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2. Relation to Existing Symplectic Integrators 

To give context to our work, we will review the development of symplectic integrators 
in the general contex of solving any evolution equation of the form 

Bin 

^ = {T + V)w, (1) 

where T and V are non-commuting operators. Classically, the evolution of a dynamical 
variable w{qi,pi) is given by the Poisson bracket (sum over repeated indices), 

d rT^_(dwdH dwdH\ 

For a Hamiltonian of the form, 

H{p,q) = ^PiPi + v{qi), (3) 
equation (2) is of the form (1) with T and V given by 

where Fj = —dv/dqi. For this case, the constituent evolution operators e^^ and e*^^ displace 
qi and Pi forward in time via 

(li^ + ePi and Pi^Pi + eFi. (5) 

In general, the evolution equation (1) can be solved iteratively via 

w{t + e)^e'^'^+^^w{t), (6) 

provided that one has a suitable approximation for the short time evolution operator e'^'-^^^^. 
If e^^ and e^^ can be exactly implemented individually, then e^^^^^^ can be decomposed to 
arbitrarily high order in the form 



Qe{T+V) 



For example, from the second order factorization 
one can construct a fourth order approximation via 

riS(6) = r(^)(?)r(^)(-.?)r(^)(?), (9) 
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where s — 2^/^ is chosen to cancel the third order error in T^^\ and e — e/{2 — s) rescales the 
sum of forward-backward-forward time steps back to e. This fourth order Forest-Ruth (FR) 
scheme (Forest and Ruth 1990) has been independently derived many times (Campostrini 
and Rossi 1990; Candy and Rozmus 1991). The above construction, while widely known 
through the work of Yoshida (1990), was first published by Creutz and Gocksch (1989). 

The middle time step in (9) is negative. This is not accidental. Sheng (1989) and 
Suzuki (1991) have independently proved that, beyond second order, any factorization of 
the form (7) must contain some negative coefficients in the set {oj, hi}. Goldman and Kaper 
(1996) later proved that any factorization of the form (7) must contain at least one negative 
coefficients for both operators. This means that for any evolution equations containing an 
irreversible operator, such as the diffusion kernel T — |V^, no algorithms of the form (7) 
is possible beyond second order. This is because (r'|e"'^^|r) oc e"^'"'"'"^^/^^"'^^ is the diffusion 
kernel only if is positive. If Oj were negative, then the kernel would be unbound and 
unnormalizable, reflecting the fact that diffusion is a time-irreversible process. Positive 
coefficients are therefore absolutely essential for solving any evolution equation having an 
irreversible component. We have shown, and will further domonstrate below, that even 
for time-reversible systems such as quantum (Chin and Chen 2002) or classical dynamics 
(Chin 1997), purely positive coefficients give rises to very efficient algorithms with excellent 
convergence properties. 

While Shcng's proof is slightly more general, Suzuki's proof provides insight into how to 
circumvent this negative cocflicient problem. The essence of Suzuki's proof is to note that, 
for example, one has the following operator respresentation of the velocity form of the Verlet 
algorithm, 

with 

H' = T + V- -e'[T, [V, T]] + -e'[V, [T, V]] + 0{e'). (11) 

In order to obtain a fourth order algorithm, one must eliminate the two double commutator 
error terms above. With purely positive coefficients ai and fej, one can eliminate cither one, 
but not both. Thus to obtain a fourth order factorization with only positive coefficients, 
one must keep one of the two double commutators. With T and V as defined by (4), the 
commutator 

simply give rises to a new force defined by the gradient of the square of the force modulus. 
This commutator can therefore be kept. The use of this force gradient is not new. Ruth 
(1983) in his pioneering paper paper has derived a third order derivative algorithm on the 
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basis of canonical transformations. However, this algorithm was universally ignored (Chan- 
nell and Neri 1996) with no follow-up discussion. In 1997 one of us (Chin 1997) noted that 
Ruth's algorithm actually corresponds to a time-asymmetric operator decomposition and 
that his force derivative is precisely the force gradient given by [V, [T, V]]. The importance 
of Suzuki's proof is that it draws a clear connection between the demands for positive time 
steps and the necessity of keeping higher order commutators, such as [V, [T, V]]. 

If one were to keep [V, [T, V"]], then there are two distinct ways of eliminating the other 
double commutator [T, [\/, T]]. Despite the lack of an operator formalism, Rowland (1991) 
also noted that the velocity form of the Vcrlct algorithm has error terms of the form (11). 
He proposed to get rid of the [T, [V, T] term by an explicit transformation and to view the 
algorithm as fourth order for a Hamiltonian with a modified potential 

V'^V+^e'[V,[T,V]]. (13) 

The rcintcrprctation of the potential is easy, the tricky part here is that since one has 
performed a transformation to eliminate [T, [y, T]], one must transform back correctly to 
preserve fourth order accurracy for the original variables. This transformation is often more 
complicated than the original algorithm. This way of eliminating [T, [V, T]] is tantamount 
to applying operators 

where e~^^^^ and e^^^^ are the initial and final transformations respectively. They can be 
further decomposed into the product form, e*^^^^ = e"^*^^e*^^. This class of "process" 
(Lopez- Marcos et al. 1996, 1997; Blanes et al. 1999) or "corrector" (Wisdom, Holman and 
Touma 1996; McLachlan 1996) symplectic algrithms has the distinct advantage that when 
one iterates the algorithm, only the kernel algorithm needs to be iterated. If there is no need 
to keep track of intermediate results, then a fourth order algorithm is possible with only a 
single evaluation of the modified force (Lopez-Marcos ct al. 1996, 1997; Wisdom, Holman 
and Touma 1996; Blanes et al. 1999). However, by its very construction, either e"*"*^^-* or e*"*^*^^ 
must contain negative time steps, and cannot be applied to equations with an irreversible 
operator. Nevertheless they are valuable addition to the symplectic repertoire for solving 
reversible dynamical problems. More recent process algorithms also uses a modified potential, 
but only in the form (13) and only in the context of a second order kernel algorithm (Blanes 
et al. 1999). Lopez-Marcos et al. (1997) have used a two-fold Rowland kernel that could 
have been fourth order, but the correct parameter values for that were never considered. 

In the operator formalism, a more direct way of ehminating [T, [V, T]\ is to add more 
operators symmetrically. In this approach, the ehmination of [T, [V, T]\ is built-in, with no 
need of an extrinsic transformation. McLachlan (1995) has done this for algorithms 4A and 
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4B (described below) in the context of a slightly perturbed Hamiltonian system. But his 
algorithms were not truly fourth order because he did not bother to keep the commutator 
[V, [T, V]]. These pseudo-higher order algorithms (McLachlan 1995; Chambers and Murison 
2000; Laskar and Robutel 2000), which only require a single error commutator to vanish at 
each order, are much simpler than regular higher order algorithms, which require all error 
commutators to vanish at each order. There is also a widely cited reference by Koseleff 
(1993) that purported to contain a fourth order scheme using modified potentials. Koseleff 's 
algorithm would have been algorithm 4A discussed below, unfortunately, his coefficient for 
the commutator [V, [T, y]] is incorrect. It should be 1/72 rather than 1/24. 



3. Fourth Order Forward Symplectic Integrators 

Suzuki, on the basis of McLachan's work (McLachlan 1995), retained the commutator 
[T, [V,T]] and wrote down two fourth order factorization schemes 4 A and 4B with purely 
positively coefficients (Suzuki 1996). He did not implement them classically or quantum 
mechanically, nor in anyway demonstrated their usefulness. One of us (Chin 1997) derived 
schemes 4A and 4B by elementary means independent of McLachan's work, and have explic- 
itly implemented them classically and demonstrated their effectiveness in solving the Kepler 
problem. Moreover, algorithms 4C (Chin 1997), which is the direct descendant of Ruth's 
third order algorithm Ruth (1983), in all cases tested, has have much smaller error coeffi- 
cients than 4A and 4B. Prior to Chin's work, we are not aware of any implementation of 
fourth order FSI for solving any problem. 

Since the commutator and the gradient force term occur frequently in the following, we 
will define 

U{t)^[V{t),[T,V{t)]] and G,(t) ^ V.(|F(t)|2) (15) 

to simplify notations. The time-dependent form of the four FSI derived in Chin and Chen 
(2002) are: 

with V defined by 

V{t)^V{t) + ^e'U{t), (17) 

corresponding to a modified force 

m^F{t) + ^e'G{t). (18) 
Given po = p(t) and qo = q(t), transcribing each operator in (16) according to (5) yields 
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the following explicit algorithm 4A for advancing the system forward to from t to t + e, 



Pi 


= Po + 


geF(qo,t) 


qi 


= qo + 


1 

-e pi 


P2 


= Pi + 


^eF(qi,i + e/2) 


q2 


= qi + 


1 

-eP2 


P3 


= P2 + 





(19) 



The last p and q are the updated variables, i.e., q(t + e) = q2 and p(t + e) = pa. For 
long time integration, when intermediate time results are not needed, the first and last force 
evaluation are the same and only need to be evaluated once. Thus algorithm 4A only requires 
two evalutations of the force and one evaluation of the force gradient to be fourth order. Note 
that forces are to be evaluated at an intermediate time equals to the sum of time steps of 
all the preceding T operators (Suzuki 1993; Chin and Chen 2002). 

Similarly algorithm 4B in operator form is 

rj4) ^^-^ ^ ^t,eT^leVit+a,e)^t,eT^^^eVit+a,e)^toeT^ (20) 

where 

= = 1(1 --L), ,. = -L, a, = l(l-±), a. = i(l + ;^). Pl) 
and with V given by 

V{t)^V{t) + Coe^U{t). (22) 

This is just a modified force with a different coefficient cq = ^(2 — y/S). The transcription 
to an explicit algorithm is similar to the 4A case and will not be repeated. Algorithm 4B 
requires two evaluations of the force and two evalutation of the force gradient. After we have 
discovered the one-parameter family of algorithms discussed below, we realize that one can 
improve 4B by eliminating one evaluation of the gradient force. This is done by concentrating 
the commutator term at the center of algorithm: 

T^f\e) = Qt'^<^T^leV(t+a2e)^ltieT^coe^U(t+e/2)^ltieT^^eV{t+aie)^toeT_ (23) 



Just to be sure that there is no confusions, we transcribe this algorithm 4B' as follows 
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qi = qo + ^oe Po 

Pi = Po+ 2^F(qi,i + aie) 
1 

q2 = qi + 2^i^Pi 

P2 = Pi + Coe3G(q2,t + e/2) (24) 
1 

q3 = q2 + 2*i^P2 

Pa = P2 + ^eF(q3,t + a2e). 
q4 = q3 + t2eP3- 

Again, the last two variables are the updated variables, q{t + e) — q4 and p{t + e) = P3. In 
contrast to algorithm 4A, 4B' only requires two force and one gradient evaluations for every 
update and can be used to produce continuous outputs. 

One of us has shown (Chin 1997) that Ruth's orginal third order algorithm (Ruth 1983) 
simply corresponds to the following asymmetric third order decomposition: 



e.(T+y)^g4Tge|yg.fTg.iy'^ (25) 

with V = V + [T, V]]. The power of the operator approach is that one immediately 

realizes that the third order error in (25) can be eliminated by symmetrization, since a 
symmetric decomposition can only have errors of even order in e. To symmetrize, take two 
third order algorithms (25) at half the time step, flip one over and concatenate. Algorithms 
4C and 4D correspond to the two possible ways of concatenation: 

r^^\e) = ^^e^^^+^^l^>)eii'^^^^+^l^e^^e^^^+^/^e^^^ (26) 

r^) (e) = gi6y(t+6)gi6Tg|.y(t+26/3)gi6Tg|ey(t+6/3)gi6Tgi6y(t) ^ ^37) 

where V{t) is the same modified force (17) used in algorithm 4A. 

Algorithms 4A, 4C and 4B' are special cases of a one-parameter family of algorithms 

T^^^,{e) = Qi3eT^V3eV{t+a3e)^t2eT^eW{t+a2e)^tieT^vieV{t+aie)^toeT (28) 

parametrized by t^. Given to, the rest of the coefficients are: 

1 11 / X / X 

^1=^2 = 2-^0, ^3 = ^0, ^1 = ^3 = g^^-32^' ^2 = 1 - (Vl + V3), (29) 
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ai — to, a2 — 1/2, as — 1 — to, with W{t) defined by 

W{t)=V2V{t)+uoe^U{t), (30) 

and ^ 

Un = — 
" 12 



1 1 

1-^ ^ + 



(31) 



l-2to 6(l-2to)3- 

All coefficients are positive for to in tlie range [0, |(1 — ^) ~ 0.21]. At to = 0, one lias 
algoritlim 4A. At to = ~ 0.167, one has algorithm 4C. At the upper limit of to = 
|(1 — ■^), V2 = 0, and one recovers algoirthm 4B'. One can therefore changes continuously 
from algorithm 4A to 4C to 4B' by varying the parameter tg. This is very useful in cases 
where any two of the above algorithms have convergence errors of opposite signs. By varying 
to, one can set that error to zero with no additional computational effort. For solving the 
gravitational three-body problem in the next section, we use to — 0.138, which is intermediate 
between that of 4 A and 4C. 

There is also a one parameter family of algorithms connecting algorithms 4B, a variant 
of 4D and 4A (Chin and Chen 2002). This family generally requires four force evaluations 
plus the force gradient. We will not consider it further in this work. 



4. The Circularly Restricted Three-Body Problem 



We compare the efficiency of these FSI by solving the planar circularly restricted 3-body 
problem defined by 

(Pr If 1 

a(r,t) = -- ai(r,t) + a2(r,t) 



where for i = 1,2, 



dt^ 



Sii{r,t) = ^ ^^^^^ and Si = \r - ri{t)\. 



(32) 



The two attractive centers ri{t) and r2{t) orbit about the origin in circles with angular 
velocity cu — 1: 

ri{t) = --(cos(t),sin(t)), 

^2it) = ^(cos(t),sin(t)). 
The gradient term is no more complicated than the force itself, 

V(|a|') = -^[Ciai + C2a2] (33) 
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with 

C2 = 2S:^^ + 3>52ai • a2 - S^^ 

We solve the problem directly in the space-fixed frame in which the force, or the accel- 
eration field as defined above, is explicitly time-dependent. The initial condition Fq = 
(0,0.0580752367) and Vq = (0.489765446,0) produce an intricate and lengthy "Chinese 
Coin" orbit, useful for testing algorithms. The actual period is 187r, but in the present 
case where the two attracting centers are indistinguishable, the orbit repeats after 97r. For 
this work we will consider the orbital period to be P = Qtt. 

In Fig. 1 we compare the trajectory after three periods {t — 277r) using the Forest- 
Ruth algorithm and the our 4B' algorithm. In Figure 2 we compare the same trajectory 
using algorithm M and algorithm 4C. M is McLachan's fourth order, four force-evaluation 
algorithm (McLachlan 1995): 

Tj^\e) = e*5f^e^^^^^*'^"*^^e**^^e''^^^^*'^"^^^e*^^^e^^^^^*"^"^^^e*2^^e''^^^^*'^"^^^e*^^^ (34) 

where Vi = 6/11, ^2 = 1/2 — vi, Oj = 5^j=i ^^id 

642 + ^471 121 



^'"^'^ 392^' 3924^^^"^^' ^3 = 1-2(^1 + ^2). 

Algorithms RF and M are representative fourth order symplectic algorithms with negative 
intermediate time steps. Both are outperformed by 4B', which only requires two force and 
one gradient evaluations. We have also tested the standard four force-evaluation Runge- 
Kutta and the three force-evaluation Runge-Kutta-Nystrom algorithm (Battin 1999). At 
the same step size, both Rung-Kutta type algorithms are unstable and cannot produce a 
bounded trajectory. Algorithm 4 A and 4B are comparable to 4B', with 4D besting even 4C 
because it uses one more force evaluation. 

The accuracy of each algorithm can be quantitatively assessed by monitoring the Jacobi 
constant. The Jacobi constant is usually given in the co-rotationg frame in which the two 
attractive centers are at rest. However, it is just as easy to transform the Jacobi constant back 
to the space-fixed frame where it has a simple interpretation. For our circularly restricted 
three-body problem defined by (32), the Jocabi constant J is given by 

J = - ^—-T - . ^— ^ - 2r X V, (35) 

|r-ri(i)| \v-V2{t)\ ^ ' 

which is twice the difference between the energy and the angular momentum. For fourth 
order algorithms, we expect J{t) — jQ-\-e^Ji{t) and the step-size independent error coefficient 
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J4{t) can be extracted from 

J,{t) ^ hm ■ 

As one decreases the step size, J4(t) converges to a fixed shape characteristic of the algo- 
rithm. In Fig. 3 we show this coefficient function for the regular fourth order Runge-Kutta 
algorithm at a step size of e = P/50000 ^ 0.00056. In the course of one period, the third 
body has five close encounters with the two attractive centers, resulting in five error spikes 
at t/P=l/10, 3/10, 5/10, 7/10, and 9/10. Between each close encounter there are also 
four minor encounters, resulting in barely discernable error blips. Fig. 3 demonstrates the 
distinction between symplectic and non-symplectic algorithm. For symplectic algorithms, 
the error recovers back to zero after each close encounter and remain bounded after each 
period. For Runge-Kutta type algorithms, each encounter produces an irrecoverable error. 
The accumulating error then grows linearly in time without limit. 

Fig. 3 shows the qualitative difference between symplectic and non-symplectic algo- 
rithms. In Fig. 4, we show the quantitative difference between forward and non-forward 
integrators by examining the detail shape of the error spike at t/P=l/10. In Fig. 4a, the 
left figure, we compare the error function J^lt) generated by algorithms FR, M and Cor with 
that of 4A. Cor is a processor or corrector algorithm of the form: 

r^:)(e)^e^(^)r,S(e)e-^(^), (36) 

where the kernel is a second order algorithm (2M) 

r2S(e)=ei^V^'(*+^/2)gier ^37^ 

with same modified force V as defined by (13). This second order kernel has been used in 
most processor or corrector algorithms (Lopez- Marcos et al. 1996, 1997; Wisdom, Holman 
and Touma 1996; Blanes et al. 1999). However, a fourth order corrector with analytical 
coefficients seems to have been overlooked in all of the above references. This corrector is 
(Chin 2002) 

gC(e) ^ QV2eV{t+lti+t2]e)^t2eT^vieV{t+tie)^tieT^ ^gg-j 
g-C(6) ^ ^-tieT^-vieV(t-t2e)^-t2eT^-V2eV{t)^ ^gg^ 



with 



2^' 2V3^' 2^3 2^/3^3' 2V3^- 

Note that at the end of e~'^^^\ time is to be set to t — {ti-\-t2)e. This corrector is unexpectedly 
short because two error terms having identical coefficients can be eliminated by a single 
corrector operator. This is the fastest fourth order algorithm when no intermediate results 
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are needed, requiring only one evaluation of the force and its gradient. When continuous 
output is needed, as it is here, it is the slowest algorithm with five force and one gradient 
evaluations. Its error function is a factor of two smaller than RF's, comparable to that of M, 
but an order of magnitude greater than that of 4A. In Fig.4b, the right figure, we compare 
all the FSI considered in this work. Algorithm 4B', with one less force gradient evaluation, 
has only half the error of 4B. The one-parameter family algorithm ACB' with to = 0.138 has 
only one-third the error of 4C. The error height of ACB' is fully two orders of magnitude 
smaller than those of non- forward symplectic integrators. In table 1, we give the inverse 
error height of each algorithm, l/\hi\, normahzed to that of RF. Thus, algorithm ACB"s 
error height is ~ 300 times smaller than that of RF and ~ 150 times smaller than that of 
algorithm M. 

To further assess the convergence properties of these algorithms, we examine the con- 
vergence of the third body's energy after one period as a function of the step size e. The 
results are plotted in Fig. 5. All can be very well fitted by single fourth order monomial Qe^ 
as shown by solid lines. One can basically read off the efficiency ranking of each algorithm 
from left to right following flatness of the convergence curve. The convergences of RK and 
RKN are similar, despite the fact that RK uses one more force evaluation. For this cal- 
culation, algorithm Cor is very efficient, requiring only one evaluation of the force and its 
gradient to best FR. But surprisingly, the second order kernel algorithm 2M seemed fourth 
order even without the use of correctors! It is nearly indistinguishable from the complete 
corrector algorithm Cor. There is a plausible explanation for this behavor. The choice of the 
modified force (13) is such that the second order error terms [T,[V,T]] and [V,[T,V]] have 
exactly the same coefficient but opposite in sign. 

picr^a/'picr ^ ^e{T+V+j^e^T,[V,T]]-^e^V,[T,V]]+-)^ (40) 

Evidently, the magnitudes of [T,[V,T]] and [V,[T,V]] arc very nearly equal and their associ- 
ated errors cancel one another. We have explicitly checked that depending on whether the 
coefficient of the modified force is less or greater than 1/24, the energy error is positive or 
negative respectively, in accordance with (40). 

As in the Jacobi constant case, all forward symplectic algorithms perform better than 
non-forward ones. Interestingly, algorithm 4ACB' is substantially better than 4A but not 
4C. Algorithm 4C and 4D appear to be anomalous; their error is definitely fourth order, but 
with exceedingly small coefficients q. In table 1, we give the inverse of error coefficient l/|ci| 
normalized to that of RF. Thus the error coefficient of RK or RKN is ten times larger than 
that of RF, while the error coefficient of 4C is nearly 2200 times smaller than that of RF. To 
give an more quantitative comparison, one should factor in the number of force and gradient 
evaluations. Although it is not unreasonable to assume that the gradient would require more 
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effort than evaluting the force, when the force and gradient are evaluted at the same time, 
the additional effort to compute the gradient, as shown by (33), is minimal. A comparison 
with RF and M, which requires three and four force evaluations respectively, should give a 
reasonable gauge of FSl's effectiveness. Taking the fourth root of l/|cj| gives the effective 
step size relative to RF for attaining the same error. Thus in the case of 4C one can use 
steps size (2200)^/*^ ~ 7 times as large as RF's, four times as large as M's and 12 times as 
large as RF's. 

Some of the observed anomalous behavior may be due to the fact that we have choosen 
an end point that is too "tame". After one period, the third body returns to a position 
that's far from both attractors with no substantial error in the Jacobi constant. To test 
this hypothesis, we compute the energy again at mid period, when the third body has a 
close encounter and the error in the Jacobi constant is at a peak. This is shown in Fig.6. 
One immediately sees that in this case the kernel algorithm 2M converges quardratically 
as it should and is distinct from the fourth order corrector algorithm Cor. Moreover, the 
convergence of both 4C and 4D is now clearly fourth order and bested by the optimized 
algorithm 4ACB'. However, fourth order fits to 4B' and 4ACB' are exceedingly unnatural, 
with very small coefficients. Both can be well fitted with a single eighth order curve as shown. 
The inverse of the error coefficient is a again given in Table 1. Despite some inexplicable 
behaviors, it is clear from Table 1 that forward symplectic integrator as a class, can be orders 
of magnitude more efficient than non-forward integrators. 



5. Conclusions 

Forward symplectic algorithms are the only composition algorithms possible for solv- 
ing evolution equations with a time- irreversible kernel (Forbert and Chin 2001; Auer et al. 
2001). However, even for reversible equations FSI have been shown to be far more efficient 
than non- forward symplectic integrators (Chin 1997; Chin and Chen 2002). In this work, 
we have further demonstrated their efficiency in solving classical dynamical problems with 
time-dependent forces. In the circularly restricted three-body problem, this class of FSI 
can have errors orders of magnitude smaller than non-symplectic or non-forward algorithms. 
While FSI can be used for any classical calculations, such as molecular dynamcis (Omelyan, 
Mryglod and Folk 2002), they are particularly ideal for doing long time, high precision evolu- 
tion of gravitational few-body problems. The force gradient is easily calculable and no more 
time consuming than evaluating the force. In comparison to other algoirithms, FSI showed 
excellent convergence behavior even at close encounters. The existence of a parametrized 
family of algorithms allows one to optimize the algorithm for individual applications. This 
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family of FSI should be applied to more realistic and more complex astrophysical problems. 

While FSI have been largely overlooked in the development of classical symplectic in- 
tegrators, they are precisely in accordance with McLachlan and Scovel's recommendation 
(McLachlan and Scovel 1996) that "derivative" symplectic algorithms should be developed. 
This work suggests that a systematic way of deriving higher order "derivative" algorithms 
is to devise factorization schemes that retain higher order commutators. 

Process or corrector algorithms can achieve fourth order accurracy with only one eva- 
lutation of the force and its gradient. If 4A or 4B' are used as kernels, one should be able 
to achieve sixth order accuaracy with only two evalutions of the force and one evaluation of 
the gradient. The use of fourth order FSI as kernel algorithms would generate a new family 
of sixth order process symplectic algorithms with minimal numbers of force evaluation. This 
work is currently in progress. If one docs not insist on having purely forward time steps, 
then many higher order algorithms can be generated on the basis of these fourth order FSI. 
See extensive constructions by Omelyan, Mryglod and Folk (2002). 

At this time, despite intense effort, no sixth order FSI has been found. At the other 
hand, no proof has been presented that this cannot be done. Clearly research in FSI is still 
in its infancy and deserves further studies. 

This work was supported, in part, by the National Science Foundation grants No. PHY- 
0100839. 

* Present address: Department of Mathematics, Texas A&M University, College Station, 
TX 77840. 
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Algorithm FR 




Algorithm 4B' 




-1 I 1 1 1 1 -1 I 1 1 1 1 

-1 -0.5 0.5 1 -1 -0.5 0.5 1 

X X 

Fig. 1. — The trajectory of the third body in the space- fixed frame after three orbits in a 
circularly restricted three-body problem. The time step used is very large, e = 97r/5000 ~ 
0.0057. FR is the Forest-Ruth algorithm which uses three force evaluations per update. 
Algorithm 4B', Eq.(23), uses two force and one force gradient evaluations. At this large step 
size, both fourth order Runge-Kutta and the Runge-Kutta-Nystrom algorithms are unstable, 
producing only chaotic trajectories shooting off to infinity. 



Algorithm M Algorithm 4C 

1 I 1 1 1 1 1 I 1 1 




.1 I 1 1 1 1 .1 I 1 1 1 1 

-1 -0.5 0.5 1 -1 -0.5 0.5 1 



Fig. 2. — Same number of orbit and time step size as in Fig.l. M is McLachan's fourth order 
algorithm which uses four force evaluations per update. Algorithm 4C, Eq.(26), uses three 
force and one force gradient. 
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Algorithm RK 



Algorithm 4C 
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Fig. 3. — The error coefficient of the Jacobi constant as computed by the fourth order 
Runge-Kutta algorithm and the forward symplectic algorithm 4C. Note the relative scale, 
algorithm 4C's error coefficient is two orders of magnitude smaller. 
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Fig. 4. — The error coefficient of the Jacobi constant at t/p=l/10 due to non-forward 
symplectic algorithms FR, M and Cor (left figure) and forward symplectic algorithms (right 
figure). Algorithm 4A is drawn with dotted lines in both for comparison. 
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Energy After One Period 




0.001 0.002 0.003 0.004 0.005 0.006 

8 

Fig. 5. — The convergence of the energy after one period. For the ease of comparion, the 
signs of 4B', RK, and RKN have been changed from negative to positive. All can be very 
well fitted by solid lines of the form Qe^. Though not visible, the convergence curves for 
4C and 4D definitely has a negative fourth order bend before turning positive. The plus 
sign corresponding to the kernel algorithm 2M, which is nearly indistinguishable from the 
complete corrector algorithm Cor denoted by hollow circles. 
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Energy at Half Period 




Fig. 6. — The convergence of the energy at mid period where the error in the Jacobi constant 
is at its peak. The sohd hnes are fourth order fits. The dotted hne is a quadratic fit to the 
kernel algorithm 2M. The dashed lines are eighth order fits to 4B' and 4ACB' whose fourth 
order coefficients are exceedingly small. The upright solid triangles are FR's results, which 
has an initial negative fourth order slope, very close to that of M, before turning positive. 
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Table 1: The inverse error height in the Jacobi constant and the inverse fourth order error 
coefficient in term of RF's value. For example, algorithm 4ACB"s maximum error in the 
Jacobi constant is 295 times smaller than that of PR's and 118 times smaller than that of 
M, McLachlan's alorithm. After one period, algorithm 4C's energy error coefficient is 2200 
times smaller than that of FR and 1100 times smaller than that of Cor. At mid period, the 
fourth order error coefficient of 4B' and 4ACB' are too small to be extracted with confidence. 
Both can be well fitted with an eighth order error term as shown in Fig. 6. 
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